function [h maxfield] = jcb_plot_near_field(filename,savefigure,fileformat)

%% initial setup
crystal_middle_size = 35;
crystal_top_size = 20;
crystal_height = 35;
num_theta_points = 100;
linewidth = 2;
draw_bowtie = false;
draw_crystal = false;
colorname = [1 1 1-1e-5];
individual_scales=false;
if nargin<3
    fileformat='png';
end
if nargin<2
    savefigure=false;
end

if nargin<1
    % get user to pick files
    [filelist,pathname] = jcb_pick_files('*nearfield*.txt', ...
    'Please select files for which to plot the near field',true,pwd);
else
    if strcmpi(filename,'all')
        filename='';
    end
    if ~exist(filename,'file')
        [filelist,pathname] = jcb_select_files(filename);
    else
        [pathname, filename, ext] = fileparts(filename);
        filelist={[filename ext]};
    end
end

%%
if ~individual_scales
    % find max field value
    maxfield=0;
    for ii=1:size(filelist,1)
        %% load file & extract data
        filename=strtrim(filelist{ii});
        profile  = load(fullfile(pathname,filename));
        maxfield = max(max(profile(:)),maxfield);
    end
end

for ii=1:length(filelist)
    filename=strtrim(filelist{ii});
    disp(['creating figure for ' filename]);
    profile = load(fullfile(pathname,filename));
    lambda=profile(1,1);
    x=profile(2:end,1);
    y=profile(1,2:end);
    offsetx=diff(x)/2;
    offsety=diff(y)/2;
    offsetx=offsetx([1:end end]);
    offsety=offsety([1:end end]);
    profile=profile(2:end,2:end);
    
    
    %% plot profile
    h = figure(ii);
    clf;
    xindices=21:length(x)-20;
    yindices=21:length(y)-20;
    jcb_contourplot(x(xindices)-offsetx(xindices),y(yindices)-offsety(yindices),log10(profile(xindices,yindices)));
    if individual_scales
        caxis([0, log10(max(profile(:)))]);
    else
        caxis([0, log10(maxfield)]);
    end
    axis off;
    colorbar off;
    %title(['near field profile at ' num2str(lambda) 'nm for ' strrep(filename,'_','\_')]);

    %% draw crystal:
    new_format=true;
    if draw_crystal
        dose=regexp(filename,'dose(?<btl>\d+)_(?<btw>\d+)_(?<gap>\d+)_','names');
        if isempty(dose)
            dose=regexp(filename,'dose(?<btl>\d+)','names');
            dose.btl=dose.btl;
            dose.btw=dose.btl;
            dose.gap=num2str(20+(80-str2double(dose.btl)));
        end
        mesh=regexp(filename,'m(?<mesh>\d+)','names');
        fileinfo.btl=str2double(dose.btl);
        fileinfo.btw=str2double(dose.btw);
        fileinfo.gap=str2double(dose.gap);
        fileinfo.mesh=str2double(mesh.mesh);
        if new_format
            loc=regexp(filename,'x(?<x>\d+p?\d*)_y(?<y>\d+p?\d*)','names');
            fileinfo.crystal_x=str2double(strrep(loc.x,'p','.'));
            fileinfo.crystal_y=str2double(strrep(loc.y,'p','.'));
            if fileinfo.crystal_x > 0
                fileinfo.crystal_angle=atan(0.5);
            else
                fileinfo.crystal_angle=pi/4;
            end
        else
            loc=regexp(filename,'loc(?<loc>\d+)sweep_d(?<d>\d+)','names');
            if isempty(loc)
                loc=regexp(filename,'loc(?<loc>\d+)','names');
            end
            fileinfo.crystal_x=0;
            fileinfo.crystal_angle=45/180*pi;
            switch str2double(loc.loc)
                case 1
                    if length(fields(loc))>1
                        fileinfo.crystal_y=str2double(loc.d);
                    else
                        fileinfo.crystal_y=10;
                    end
                case 2
                    fileinfo.crystal_y=20;
                case 3
                    fileinfo.crystal_y=40;
                case 4
                    fileinfo.crystal_y=108;
                case 5
                    fileinfo.crystal_x=75;
                    if length(fields(loc))>1
                        fileinfo.crystal_y=str2double(loc.d);
                        fileinfo.crystal_angle=atan(0.5);
                        if fileinfo.crystal_y==56
                            fileinfo.crystal_y = 58 -4*sin(fileinfo.crystal_angle);
                            fileinfo.crystal_x = fileinfo.crystal_x -4*cos(fileinfo.crystal_angle);
                        elseif fileinfo.crystal_y==54
                            fileinfo.crystal_y = 58 -8*sin(fileinfo.crystal_angle);
                            fileinfo.crystal_x = fileinfo.crystal_x -8*cos(fileinfo.crystal_angle);
                        end
                    else
                        fileinfo.crystal_y=71;
                    end
                case 6
                    fileinfo.crystal_y=[];  % as no crystal used
            end
        end




         % check profile height
        fid=fopen(fullfile(pathname,filename),'r');
        fileinfo.height=fscanf(fid,'%c',300);
        fileinfo.height=regexp(fileinfo.height,'height z=(\d+\.?\d*)','tokens');
        fileinfo.height=str2double(fileinfo.height{1});
        fclose(fid);

        % calculate the crystal cross-section size at the given height
        cryst_size= crystal_middle_size ...
                    - abs(2*fileinfo.height/crystal_height - 1)...
                      * (crystal_middle_size-crystal_top_size);

        % define crystal coordinates, centred at origin
        coords = 0.5*[-cryst_size, -cryst_size;
                  -cryst_size,  cryst_size;
                   cryst_size,  cryst_size;
                   cryst_size, -cryst_size;
                  -cryst_size, -cryst_size]; 

        % rotate crystal
        for kk=1:size(coords,1)
            coords(kk,:) = ([cos(fileinfo.crystal_angle), -sin(fileinfo.crystal_angle);
                            sin(fileinfo.crystal_angle),  cos(fileinfo.crystal_angle)]*coords(kk,:)')';
        end

        if ~isempty(fileinfo.crystal_y)
            % shift coords to the correct position for the crystal
            coords = coords + ones(size(coords,1),1)*[fileinfo.crystal_x fileinfo.crystal_y];
            line( coords(:,1), coords(:,2),'color',colorname,'LineWidth',linewidth);
        end
    end

    %% draw bowtie
    if draw_bowtie
        p  = 2*atan(fileinfo.btw/2/fileinfo.btl);
        q  = (pi-p)/2;
        Rc = 9;
        L  = fileinfo.btl;
        W  = fileinfo.btw/2;
        fileinfo.gap=fileinfo.gap;
        L=L+Rc*(1/sin(p/2)-1);
        W=W+Rc*(1/tan(q/2)-1);

        inset=Rc/sin(p/2);
        Wa = Rc*cos(p/2);
        Ha = Rc*sin(p/2);

        e  = Rc/tan(q/2);
        f  = e*cos(q);
        g  = e*sin(q);

        lineX = [ [inset-Ha, L-g];
                  [L, L]
                  [L-g, inset-Ha]] - (inset-Rc)+fileinfo.gap/2;
        lineY = [[Wa, W-f];
                 [W-e, -W+e]
                -[W-f, Wa]];

        circleX = zeros(3,num_theta_points);
        circleY = zeros(3,num_theta_points);
        theta=linspace(2*pi-(pi/2+p/2),pi/2+p/2,num_theta_points);
        circleX(1,:) = Rc*cos(theta)+Rc+fileinfo.gap/2;
        circleY(1,:) = Rc*sin(theta);
        theta=linspace(pi/2+p/2,0,num_theta_points);
        circleX(2,:) = Rc*cos(theta)+L-inset+fileinfo.gap/2;
        circleY(2,:) = Rc*sin(theta)+W-e;
        theta=linspace(0,-pi/2-p/2,num_theta_points);
        circleX(3,:) = Rc*cos(theta)+L-inset+fileinfo.gap/2;
        circleY(3,:) = Rc*sin(theta)-W+e;
        fullX=[circleX(1,:) lineX(1,:),circleX(2,:) lineX(2,:),circleX(3,:) lineX(3,:)];
        fullY=[circleY(1,:) lineY(1,:),circleY(2,:) lineY(2,:),circleY(3,:) lineY(3,:)];           
        line(fullX,fullY,'Color',colorname,'LineWidth',linewidth);
        line(-fullX,fullY,'Color',colorname,'LineWidth',linewidth);
    end

    %% save figure to file
    if savefigure
        print(['-d' fileformat],regexprep(filename,'.txt',''));
    end
end

if nargout < 1
    clear h
end